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PROGRAMS FOR HIGH-SPEED FOURIER, MELLIN 
AND FOURIER-BESSEL TRANSFORMS 


D, K. Tkhabisimov, A. S. Debabov, B. I. Kolosov, 
D. A. Usikov 


We present a description of program modules for 
performing one-dimensional and two-dimensional discrete 
Fourier transforms, Mellin and Fourier-Bessel trans- 
forms, and also programs for realizing the algebra of 
high-speed Fourier transforms on a computer. The pro- 
grams developed can be used to perform numerical harmonic 
analysis of functions, to synthesize complex optical fil- 
ters on a computer, modelling the holographic methods of 
processing images. 

The programs are written in FORTRAN and are in- 
cluded in the library of programs of the video images 
processing unit "SOFI Kl in the Institute of Space 
Sciences of the USSR Academy of Sciences . 


INTRODUCTION 


/ 3* 


The need for numerical harmonic analysis arises when solving integral 
equations of convolution type. Here all the functions are given in discrete 
form. Let fm be the file of the readings for the initial function fez}' 
which is nonzero on the interval (aJJ of the real axis. A pair of discrete 
Fourier transforms (DFT) of the file f ( W has the form: 




and 


( 1 ) 


_ 1 + ft/-i JC 

{in) = D nm " - f(m . ), 


( 2 ) 


Numbers in the margin indicate pagination of original foreign text. 


1 



where N is the dimension of the files f (n) and f (m) . By means of the trans- 
forms (1) and (2) we can compute continuous Fourier transforms and the coeffi- 
cients of the Fourier series of the initial function [ 1] . The Mellin transform 
of the function is defined by the formula 12]: 

f M ffeO ~ J 31 ^ If C*D Jx . C3) 

On making the substitution x-e? under the integral sign in (3) , we obtain 

^ [Wef C4) 

where is the operator of the continuous Fourier transform. From 04) it 

is evident that the Mellin transform is carried out with the help of Cl) • The 
Fourier-Bessel function -f.fa) has the form: 


ts x r/fle -f(x) *= f £ /fx)x c/x . C5) 

In II] a procedure is shown for computing (5) with, the help of the transforms 
Cl)-C3) . Applying the Mellin transform C3) to C5) , we obtain 
y(h>0) 

/rtu/,* /k w = M 0fJ> { L (p) 


/4 



C6) 


C7) 


a . c 

The coefficient //?) is computed at the points where T ~ 

The case n = O' reduces to the case n = 1 if we note that 

_ _ *0 ' ' 

1 1 /<?>/> Jf> = t 1 T(f) 1 Jp.. 

where 
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r<7» = J- / 1 


- f j f($) S'c/$, 


To calculate (7) we can use the representation I2J; 


where t'iji '/ . When we make use of the recursion formula: 

nztiizz rM. 


Bessel's functions are computed hy means of Cl ) > since 


18) 


1(f)* 


2?r t ‘ 

fr l C £ && 


C9) 


In the present paper, programs are described which realize the transforms 
(1) , (2) , (3) and C5) . A description is also given of programs for the algebra 
of high-speed Fourier transforms, programs for calculating the Bessel and gamma 
functions, etc. The text of the programs appears in the appendix. 


1. Program for the High-Speed Fourier Transform 
1.1 The One-Dimensional Version 

The programs DQFT1 and KQFT1 are designed to carry out direct and inverse 
DFT of one-dimensional files described in the user's program as a complex. In 
accordance with the principles of organization of the complex "SOFI" I4J, access 
to the programs is accomplished as follows : 

1) In the basic program the user defines 

1 tmu EX" JteJ 

CO two V F2CT (3 M, 11/3), d 

>fi/=zfv 

. • I 

Illegible in original foreign text. 
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where A is a working one— dimensional file with, the DFT program to he processed 

2) The file to be processed is entered into the working file A and accesses 
one of the programs of the high-speed Fourier transforms, e.g. , CALL DQFT1; 

3) The Fourier transform of the initial file is placed in 'the file A; the 
error in the calculation does not exceed 10 when n&Q. 


1.2 The Two-Dimensional Version 


The direct and inverse DFT of two-dimensional files is performed with the 
help of the programs DQFT2 and RQFT2. Here not every variable undergoes (1) and 
(2) . In the basic program it is necessary to write: 


Wn, fj cross) J Tg(s)i ^ B 


/6 


where A, B are the one- dimensional working files with the DFT program to be 
processed ,{//.= J2 , il a ~2 Moreover, analogous to Sec. 1.1, a description of 
the programs used in the high-speed Fourier transform programs can be found in 
13]. 


When formulating the basic program, the user must designate the library 
of the complex for which the card 

//,, PAEEE tJ X 

must be inserted after the card 

7/~ T0& _ A’ a WE 


2. The Fourier-Bessel Transform 

The Fourier-Bessel coefficients C5) are computed by the program 

f — 

FTESjEj VTJ-A/8 a WJpjTM)- Here F is the file for the readings of the 
initial function f(x) at the points ^ ~ , NT is the number of readings on 

the interval *o5 WE is the subscript of the Bessel function, W is the 
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file of the function • 7/5^ v TM = t • For access to the program FTFB, 

• ~ max 

the user must write in the basic program 


Co'MP/rx 

HQ/vL t icjjj 1 3J, A _ 


where m is the dimension of the file. The result is located in the file F. 

The values of the function W. (w) are calculated by means of the program 
FT/ T (v/j // /) T/Sj> T 7-f . Here W is the desired file of values of the function /7 

is the mnnber of values, M , F /- ^ ~ ^ , the function 
is computed at the points ; F.- O f . y / TTff/~ . In the program FILT, the 

gamma function is calculated with the help of series (cf. Sec. 3); however, it 
is also possible to make use of the representation (8) * The programs FTFB and 
FILT are accessed by means of the operator CALL. 


3. Calculating the Bessel and Gamma Functions 

The program BESFCBF,TF,NF) is calculated by means of formula C9) NF of the 
Bessel functions at the point TF: 

X (TF), j; (TF), ■ ■ v / f -, (TF) 

and the result is placed in the file BFCNF) . 

The program C/F //f } T0 t J}j, computes the function at 

the points where (f/j- and F0 are the upper and 
lower integration limits in (8) . The result is located in the file G(NT) . Be- 
fore accessing programs BESF and GAMMA1, it is necessary to describe: 

Co np / ex T of 

con no v Pict?35?) ; ^ 


The program GAMMA (Z) calculates the gamma function at the point // 
by means of the formula [5] : 

T (z) = mT / 7 -^ 4 ~ ■ 


K. 


£=/ 
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The result is located in Z. 


4. Auxiliary Programs 

Below we give descriptions of programs whose use simplifies access to the 
high-speed Fourier transform programs and accelerates the process of formulating 
a basic program. All programs work with complex files. 

(j?j Sj K , /. )- the program selects line number K from file 

AC_NX,NY) and when L = 0, copies the information from file A into file B(NY); when 
L = If it copies the information from file B(NY) into file A. 


2) C./.M A/ (dj 3j A/X / JCj Aj " the program selects column number K from 

file A(NX,NY) ; when L = 0, it copies the information from file A into file B(NX); 
when L - 1, it copies the information from file B(NX) into file A. 


3) Q U T C OAjfjQj //X ; /V^ — the program copies the information from the 

COMMON-hlock of element K in file C into file ACNX^NY) * 

4) TO.COM Ms — the p r0 g ram copies the information into the 

COMMQN-block of element K in file C from file A(NX ? NY) • 

5) TMU/TfT&j 0 ; MX,A/yJ~ the program multiplies out the two two-dimen- 
sional files A(NX ? NY and B(NX ? NY), component by component. The result is located 
in the file SCNX,NY). 


[6)1 TH//yp(Jj & //Xj A/yJ " the program transposes 
The result is located in the file B(NX,NY). 


the matrix A (NX, NY) . 


] =T)1 . COA/T G (Aj Bj x/yy — the program forms the complex conjugate of 

the matrix A (NX, NY) . The result is located in the file B(NX,NY), 

mo vi & (a, vx hv, a/Xj WjB; A/iA/pw, ^this is a program for the linear in- 
terpolation of a function, given in Cartesian coordinates by the file of values 
A(NX,NY) on a network with sampling steps HX and HY along the axes OX and 0Y, 
respectively, into the polar exponential network (.£ Jp 5 '") > where 

-cdsiAx.00, 

Here T" P/rf **”'**-, /At -’ — number of subdivisions of the network with respect _/9_ 
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to the variable t, NF — number of subdivisions of the network with respect to 
the variable fit. The result is located in the file B (NT,NF) . 

9) ZEVolfctyMj/t A/X, A/y t }l /;J — is a program of linear interpolation 

of the values of a function given in the Cartesian file iF^A^A/i) into the 
Cartesian network rotated through the angle AL relative to the origin of coordi- 
nates. The origin of coordinates corresponds to the point -r d 

NX and NY are the numbers of the nodes of the new network with respect to each 
variable (the sampling step is not changed); V(NX,NY) is the file of values of the 
initial function in the new network. 

10) COR. F ( Fj /y'y± - the program calculates the correlation 

function of the too images in the plane given by the files F(NX,NY) and V(NX,NY) 
of the form: 

= X4 c/sc , ao) 

x<rhere 3? . ^ are vectors with components and respectively. In 

the result, the file S(NX,NY) of the readings of the correlation function is ob- 
tained . 


, CONCLUSION 


In the present paper the programs described are included as object modules 
in the complex '‘SOFI" anil are intended for various problems in spectral processing; 
recognition of images, refinement of images, the photographing of smog, the effects 
of blurring, for various statistical and physical problems required when applying 
the algorithms for the Fourier, Mellin and Fourier-Bessel transforms. 

The authors express their deep appreciation to V. G. Zolotykhin for useful 
discussions . 
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